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ABSTRACT 


In  this  thesis,  the  use  of  the  Wentzel-Kramers-Brillouin  (WKB)  Theory’  to  obtain  the 
solution  to  the  Helmholtz  Equation  governing  the  acoustic  normal  modes  is  examined. 
Specifically,  uniformly  valid  WKB  solutions  for  four  classes  of  acoustic  normal  modes 
in  the  ocean  are  derived  and  the  accuracy  of  the  WKB  approximation  is  tested  against 
some  exact  solutions.  It  is  found  that  this  inherently  high  frequency  technique  has  an 
appreciable  accuracy  even  at  a  frequency  of  1  Hz.  A  product  of  this  thesis  is  a  computer 
program  that  solves  for  the  WKB  modes  for  an  arbitrary  sound  speed  profile. 
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THESIS  DISCLAIMER 


The  reader  is  cautioned  that  computer  programs  developed  in  this  research  may  not 
have  been  exercised  for  all  cases  of  interest.  While  every  effort  has  been  made,  within 
the  time  available,  to  ensure  that  the  programs  are  free  of  computational  and  logic  er¬ 
rors,  they  cannot  be  considered  validated.  Any  application  of  these  programs  without 
additional  verification  is  at  the  risk  of  the  user. 
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I.  INTRODUCTION 


A.  BACKGROUND 

To  model  sound  propagation  in  the  ocean  we  can  use  several  theories,  namely  Ray 
Theor}',  Parabolic  Equation  Approximation  and  Normal  Modes.  The  Ray  Theory  sol¬ 
ution  is  an  asymptotic  geometric-optics  solution,  obtainable  using  simple  ray  tracing 
techniques.  It  provides  a  simple  physical  description  on  how  sound  is  transmitted 
underwater.  However,  it  neglects  sound  diflraction  and  thus  needs  corrections  near 
caustics  and  turning  points.  Such  corrections  can  be  very  complicated  mathematically. 
The  Parabolic  Equation  method  is  less  physical  but  is  a  full-wave  solution.  An  asset  of 
Parabolic  Equation  is  its  capability  to  handle  variable  bottom  bathymetry  well.  Its  lim¬ 
itation  is  that  it  does  not  accurately  model  sound  energy  propagating  in  steep  angles. 
Normal  Mode  Theory  gives  a  physical  full-wave  solution  to  the  wave  equation.  It  has 
some  computational  difficulties  in  handling  range  varying  sound  speed  fields. 

A  computational  difficulty  associated  with  normal  mode  models  applied  to  a  range- 
dependent  ocean  is  that  the  normal  modes  associated  with  a  great  number  of  profiles 
must  be  computed  and  the  results  stored.  Large  storage  and  processing  time  are  re¬ 
quired  if  straight-forward  numerical  methods  such  as  finite  differences  are  used  to  com¬ 
pute  the  normal  modes  (Chiu  and  Ehret,  1990).  The  use  of  finite  differences  methods 
requires  the  discretized  mode  functions  and  their  eigenvalues  at  a  great  number  of  points 
to  be  stored.  Moreover,  at  high  frequencies  these  methods  require  the  computation  of 
eigenvectors  and  eigenvalues  of  large  matrices,  which  can  result  in  significant  increases 
of  processing  time  and  numerical  noise  in  the  solution.  In  this  thesis  we  examine  an 
asymptotic  expansion  method  that  has  the  potential  to  overcome  this  dilTiculty.  The 
method  is  called  Wentzel-Kramers-Brillouin  (WKB)  theory. 

B.  THE  NORMAL  MODE  APPROACH 

The  wave  equation  governing  the  sound  pressure  p  in  the  ocean  is 


where  c  =  c(z;  r,  6)  is  the  sound  speed  and  z  is  the  vertical  coordinate,  r  the  range  and  9 
the  azimuthal  angle.  In  the  three-dimensional  coupled  mode  model  of  Chiu  and  Ehret 
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(1990),  the  pressure  is  expressed  as  a  linear  combination  of  the  local  normal  modes  Z„ 
such  that 


oo 

pir,  d,  z)  =  '^Z„iz;  r,  d)P„{r,  e)T(t)  . 

For  a  harmonic  frequency  time  dependence 

T{t)  =  c'"' 


where  cu  is  the  acoustic  angular  frequency,  the  local  modes  at  each  horizontal  location 
(r,  6)  are  required  to  satisfy 


dz^ 


=  0 


and  the  appropriate  boundary  conditions.  This  equation  is  usually  known  as  the 
Helmholtz  Equation.  The  constant  k,  is  the  horizontal  component  of  the  wavenumber 
vector  whose  magnitude  is  given  by 


In  general,  there  are  many  possible  values  of  k„  (eigenvalues)  satisfying  the  Helmlioltz 
Equation.  For  each  k„  there  is  an  associated  mode  Z„  (eigenfunction).  In  a  range- 
dependent  sound  speed  field  this  eigenvalue-eigenfunction  problem  must  be  solved  at  a 
great  number  of  horizontal  grid  points. 

C.  THESIS  OBJECTIVES  AND  OUTLINE 

The  main  objective  of  this  thesis  is  to  examine  the  use  of  the  WKB  theory  to  solve 
the  Helmholtz  Equation  governing  the  acoustic  normal  modes.  This  includes: 

1.  the  development  of  the  various  WKB  formulae  for  the  four  classes  of  normal  modes 
which  can  exist  in  single  duct;  channel  environments, 

2.  the  parameterization  of  each  class  of  normal  modes  using  a  minimum  number  of 
parameters, 

3.  the  quantification  of  errors  in  the  WKB  solution  through  comparisons  to  three 
exact  solutions. 

A  product  coming  out  of  this  thesis  is  a  computer  program  solving  for  the  mode 
parameters  for  an  arbitrary  sound  speed  profile.  The  incorporation  of  this  code  in  the 
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three-dimensional  coupled  mode  model  of  Chiu  and  Ehret  (1990)  is  expected  to  result  in 
significant  savings  of  processing  time  and  computer  storage. 

In  Chapter  II,  the  WKB  formulae  are  developed.  The  four  types  of  normal  modes 
are  discussed.  In  Chapter  III,  the  normal  modes  and  their  eigenvalues  for  three  abstract 
sound  speed  profiles  are  solved  exactly.  The  WKB  results  are  then  compared  to  the  exact 
solutions  for  an  error  analysis.  Conclusions  are  given  in  Chapter  IV. 
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II.  WKB  NORMAL  MODES 


A.  FIRST  ORDER  WKB  APPROXIMATION 
1.  The  Mathematical  Problem 

The  WKB  approximation  to  the  solution  of  a  differential  equation  is  an 
asymptotic  expansion  method  applicable  to  a  large  range  of  equations.  It  is  named  after 
Wentzel,  Kramers  and  Brillouin  who  used  it  separately  but  at  about  the  same  time  in 
1926.  However  the  principle  of  this  technique  was  developed  by  Liouville  and  Green  in 
1837.  It  was  also  used  by  Rayleigh  (1912),  Gans  (1915)  and  Jeffreys  (1924),  among  oth¬ 
ers.  The  WKB  method  is  also  known  as  the  Liouville-Green  or  WKBJ  approximation 
(the  letter  J  is  used  to  honour  Jeffreys'  contribution). 

To  compute  the  acoustic  normal  modes  the  equation  to  be  solved  is  the 
Helmholtz  Equation 


+•<12.  =  0 

dz 


(1) 


where 


C 

is  the  vertical  wavenumber,  Z,  is  the  normal  mode,  z  is  the  vertical  coordinate,  w  the 
acoustic  angular  frequency,  c=c(z)  the  sound  speed  and  k,  the  horizontal  wavenumber. 
With  a  pressure  release  surface  and  a  rigid  bottom  the  appropriate  boundary  conditions 
are 


Z„  =  0  at  the  surface  (2) 

dZ„  ... 

— - —  —  0  at  the  bottom  .  (3) 

dz 

Normal  modes  subjected  to  general  boundary'  conditions  are  discussed  in  Appendix  B. 
Note  that,  in  (1)  through  (3),  we  have  supressed  the  dependence  in  range  which  has  no 
effect  on  the  local  modes.  Equation  (1),  together  with  (2)  and  (3),  is  a  Sturm-Liouville 
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problem  where  the  Z„'s  are  the  eigenfunctions  and  the  k:„'s  are  the  eigenvalues.  The 
eigenfunctions  can  be  normalized  using  a  normalization  constant  C„  such  that 

=  QZ„  (4) 

where 


Note  that  the  integration  is  over  the  entire  depth  in  (5),  is  the  Kronecker  Delta,  and 
Z„  and  Z,  are  orthogonal  for  «  m. 

There  are  several  ways  to  obtain  the  first  order  WKB  solution.  A  physical  ap¬ 
proach  is  to  consider  ;lie  transmission  of  plane  waves  through  a  layered  medium.  The 
WKB  solution  is  obtained  as  we  let  the  the  layer  thickness  approach  zero.  This  physical 
approach  will  be  discussed  next. 

2.  Physical  Approach 

In  each  isospeed  layer  the  solution  to  the  Helmholtz  Equation  is  given  by 

z{  oc 

where  is  the  vertical  w'avenumber  in  the  j’’'  layer.  The  transmission  coefficient  from 
layer  j  to  layer  j  +  1  is  given  by  (Kinsler  et  al.,  19S2) 


T  —■} _ _ 


^yV+i  ~ ' 


i+T— 


where  Ak:^,  =  With  — 1  we  get 


^jj+l  -  ^  ^  'fin 


By  letting  1  Zi  I  =1  vve  have 
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where  is  the  value  at  the  interface  of  the  layer  and  the  (/  +  1)'^  layer  and  Az^  is  the 
thickness  of  the  m'*'  layer.  Taking  the  thickness  of  the  layers  to  zero  we  have 


or 


V  '^2n 


where  k,  is  the  vertical  wavenumber  at  z  =  z.  and 


<t> 


=  J  \K,„\dz  . 


The  vertical  wavenumber  k„  generally  can  take  on  real  or  imaginary  values.  If 
K„  is  real,  i.e.,  kI,  >  0,  the  solution  is 


^  ,  a'e'Ub’e-‘^ 

2„(z)  - - = - 

V  ^2n 


or 


a  cos  <^+h  sin  (f> 


(6) 


where  a'  and  b'  or  a  and  b  are  constants  to  be  determined  by  normalization  and  the 
boundary  conditions.  On  the  other  hand  if  k„  is  purely  imaginary,  i.e.,  <  0,  we  have 


Zn{2)  = 


ae‘^+be'~'^ 

V'RJ 


(7) 
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These  last  two  formulae  are  precisely  the  VVKB  solution  to  the  Helmholtz  Equation  in 
regions  not  close  to  a  point  where  k„  =  0  (i.e.,  a  turning  point). 

3.  Basic  Formulae 

By  substituting  the  WKB  solutions  (6)  or  (7)  in  the  Helmholtz  Equation,  one 
can  easily  show  that  they  are  exact  solutions  if 


r(r)  = 


3  1 


dK,„ 


I  1 


d^K, 


=  0 


Therefore,  if  |r(z)|<|I  the  WKB  solutions  (6)  and  (7)  are  good  approximations  to  the 
exact  solution  for  kI„  >  0  and  k^,  <  0,  respectively.  In  Appendix  A  a  more  detailed  and 
mathematical  derivation  of  the  WKB  solution  and  its  validity  is  pre.sented.  Let  us  call 
these  solutions  the  first  order  WKB  approximation.  But  (6)  and  (7)  are  not  valid  when 
=  0  (i.  e.,  at  a  turning  point)  or  even  in  a  region  where  k]„  ~  0.  Wc  need  a  solution 
valid  at  and  near  the  turning  point  (i.e.,  in  the  critical  region)  to  make  the  liaison  be¬ 
tween  the  oscillatory  region  (k:^„^0)  and  the  exponential  region  (k:^„<^0). 

Let  us  suppose  that  there  is  a  turning  point  at  z  =  0  and  k^„>0  for  z  >0  and 
<  0  for  z  <  0.  Consistent  with  a  first  order  approximation,  let  us  assume  that  near  the 
turning  point  k^„  =  yz  where 


y  = 


With  a  change  of  coordinate 


r  >/3 

C  =  -y  z  , 


the  Helmholtz  Equation  becomes 


d:^ 


-CZ„  =  o 


which  is  the  Airy  Equation  with  the  general  solution  given  by 

Z,(C)  =  . 

The  Airy  Functions  can  be  expressed  as,  with  C  >  0  , 
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..(;)=4['-./3(t=’'0-''»(t<"0] 

ara = Vt  [/-., .(I  {’'=)+/, 


where  J's  are  Bessel  Functions  and  I's  are  Modified  Bessel  Functions. 

Now  we  have  a  complet  set  of  first  order  approximate  solutions  covering  the 
entire  range  of  kJ„.  Summarizing,  we  have: 

(a)  in  the  oscillatory-  region  (k^^O) 


a,  sin  4>+hi  cos  d> 

Z„{z)  - - = - 

\ 


(8) 


(b)  in  the  critical  region  (k^  ~  0) 

K-n  =  yz 


r 

C  =  -y  '  7 


Z„(7)  =  a^AiiQ+b^DiiO  , 


(c)  in  the  exponential  region  {k]„4.0) 


Zn(^) 


a^e  ’^+bje'^ 
V'lK^nl 


(10) 


To  assure  continuity  in  Z^z)  at  the  boundaries  between  regions,  we  must  relate 
the  constants  a,  and  fe,  carefully.  An  asymptotic  expansion  of  the  critical  region  solution 
for  large  positive  values  of  kJ,  (or  — C  >  1)  must  match  the  solution  in  the  oscillatory  re¬ 
gion  and  for  large  negative  values  of  (or  C  >  1)  must  match  the  solution  in  the  expo¬ 
nential  region.  The  connection  formulae  between  the  three  regions  are  given  by 
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*- aAiXO  -> — =L=r- sm((^+ ) 


—  ^  <-  aBi(C)  — ;==-  cos(<^+  ) 

\  I  I  V  ^r/1 

where  the  left  sides  match  the  Airy  Functions  for  Kj„<^0  and  the  right  sides  match  for 
(Bender  and  Orszag,  1978)  with 

Therefore,  (8),  (9)  and  (10)  can  be  recast  respectively  as 

(a)  oscillator^’  region  solution 

a  sin(d)+  -j )+/?  cos(4>+  -j ) 

Zniz)  = - ^  ,  (ll.«) 


(b)  critical  region  solution 

Z,(2}  =  aaAiaHtcrBKQ  , 

(c)  exponential  region  solution 


Zni^) 


{alDe-^+be^^ 

V  I  Xzn  I 


ill.b) 


(ll.c) 


The  constants  k„,  a  and  b  are  determined  by  normalization  and  the  boundary 
conditions.  The  equations  for  the  eigenvalues  k„,  i.e.,  the  characteristic  equations,  in¬ 
volve  the  solution  in  the  oscillatory  region  (U.a),  as  will  be  discussed  later. 

Equations  (ll.a,b,c)  are  not  easy  to  use  in  the  computation  of  normal  mode 
shapes.  Where  does  one  stop  with  one  formula  and  start  with  another?  This  problem  is 
avoided  if,  after  the  determination  of  the  eigenvalue  k„  and  the  constants  a  and  b,  the 
computation  of  the  normal  mode  is  done  using  the  Langer  Formula  (Nayfeh,  1973; 
Bender  and  Orszag,  1978) 


Zn{z)  =  2^fn 


1 


r 


aAi 


(12) 
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with 


and 


s  =  -2l\z\  . 

This  formula  has  the  advantage  of  giving  a  single  and  continuous  solution  for  the  entire 
range  of  Bender  and  Orszag  (197S)  have  shown  that  for  and  K^  =  yz  the 

formula  gives  exactly  (ll.b),  for  Ki,^0  it  asymptotically  approaches  (ll.a)  and  for 
it  approaches  (1  l.c). 

In  rectrospect,  in  order  to  obtain  the  first  order  WKB  solution,  an  algorithm 
must  include  the  following  steps: 

(a)  application  of  a  boundarv'  condition  to  (ll.a,b,c)  to  get  the  relationship 
between  the  constants  a  and 

(b)  application  of  the  other  boundai7  condition  to  (ll.a)  to  get  the  character¬ 
istic  equation; 

(c)  with  (12)  compute  Z„(z); 

(d)  with  (4)  and  (5)  normalize  Z„(r). 

4.  Comments 

The  WKB  solution  is  a  high  frequency  approximation.  This  means  that  it  gets 
better  as  the  frequency  gets  higher.  It  must  be  noted  that,  although  Normal  Modes  are 
a  full-wave  exact  solution  to  the  wave  equation,  the  WKB  modes  arc  approximate  sol¬ 
utions  and  their  accuracy  is  frequency  dependent. 

B.  DETERMINATION  OF  WKB  MODE  PARAMETERS 

Now  that  we  have  a  uniformly-valid  first  order  solution  to  the  Helmholtz  Equation, 
we  are  ready  to  apply  the  boundary  conditions  to  get  expressions  for  the  k„'s  and  the 
constants  a  and  b  (or  their  equivalents).  There  are  four  classes  of  normal  modes  to 
consider.  Each  class  has  different  mathematical  expressions  for  the  mode  parameters  (a, 
b  ,  or  their  equivalents,  and  k„).  Therefore,  for  an  arbitrary  sound  speed  profile  the  first 
procedure  for  normal  mode  calculation  is  to  determine  which  class  each  mode  falls  into 
and  then  go  through  the  steps  described  at  the  end  of  the  previous  section. 
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I.  General  WKB  Formulae 

Wc  start  by  deriving  the  more  general  formulae  which  are  applicable  to  modes 
that  do  not  have  turning  points  in  close  vicinity  of  the  boundaries. 
a.  Class  I:  Pressure  Release  -  Turning  Point 

First  let  us  consider  a  mode  whose  oscillatory  region  is  bounded  by  a  turn¬ 
ing  point  at  the  depth  z  =  i  and  the  surface  (z  =  0).  With  the  bottom  at  z  =  //,  the 
boundary  conditions  are 

Z„(0)  =  0  (13) 


dZ„ 


(14) 


We  will  call  this  class  of  modes  as  PR-1  P  (Pressure  Release  -  Turning  Point). 

In  the  exponential  region  (2<^z  <  H),  the  WRB  solution  can  be  recast  as 


Z„  oc 


coshC0-<^^) 


V 


(15) 


with 


In  order  to  satisfy  the  boundary  condition  (14).  0^  must  be  given  by 


We  need  to  connect  (15)  to  the  solutions  in  the  critical  and  oscillatory  regions  next.  After 
application  of  the  connection  formulae,  the  results  in  the  three  regions  are: 

(a)  exponential  region 

Z„  =  — _L=-  cosh(0-0f,)  ,  (16.a) 

V  I  Kzn  I 


(b)  critical  region 


=  y(^-z) 
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Z„  =  ^‘^‘0/liK)+-^cT5i(O  . 


(I6.fc) 


(c)  oscillatory  region 


V  ^2n  ^2n 


-flK, 


in  this  region.  The  corresponding  Langer  Formula  that  asymptotically  matches 
(16.a,b,c)  is 

Z„  -  27^ (  f  I  *  I )  ■'«  {/•.!,[.(  I  «)=«]+  41  b[s{  I  ^ )>'’]}(I7) 


As  the  oscillatory  region  solution  must  satisfy  the  surface  boundary  condi¬ 
tion,  we  have 


T  )+  f  ]  =  0  . 
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K,„dz  =  (^n- —  ^n-  tan 

This  is  the  characteristic  equation  for  PR  -  TP  normal  modes.  The  solutions  to  the 
characteristic  equation  give  the  eigenvalues  k,'s. 

b.  Class  II:  Turning  Point  -  Rigid  Bottom 

The  next  class  of  normal  modes  to  be  considered  has  the  oscillatory’  region 
between  a  turning  point  and  the  bottom.  It  will  be  called  TP-RB  (Turning  Point  -  Rigid 
Bottom).  For  this  class,  let  us  express  the  solution  in  the  exponential  region, 
as 


oc  —  sinh((/)— 

sj  I  I 


with 


and 


4> 


=  f  \’<2n 


\dz 


A 


,\dz. 


U  sing  the  connection  formulae,  we  get  for  the  three  regions: 
(a)  Exponential  Region 


Z„  =  - 


1 


V  I  I 


sinh(0— 


(b)  Critical  Region 


=  yi2-z) 


C  =  -y'%-z) 
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(c)  Oscillatoi7  Region 


V  ^2/1 


(19J 


with 


in  this  region  (2<^z  <  H).  The  appropriate  Langer  Formula  for  this  class  is 


where 


A 

z—z 

1  2-Z  I 


and 


<t> 


-I' I"- 


\dz. 


Equation  (19)  must  satisfy  the  bottom  boundary  condition  (14),  therefore,  we  have 


sm 


_  2 


4*  tan 


/J 


=  0 


(21) 


where 


D  = 


1  \ 
2kL  dz 


(22) 


It  follows  from  (21)  that  the  characteristic  equation  for  TP  -  RB  is 
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for  D  <  0,  and 


k,„^/2  =  ^rt-j^Tt-tan  ’ 


<?~'^‘  _  a 
— r - Dc  ' 


for  D  >  0. 

c.  Class  III:  Turning  Point  -  Turning  Point 

The  third  class  of  modes  has  the  oscillatory  region  between  two  turning 
points,  at  2  =  2,  and  2  =  Zj  with  ij  <  z^.  They  will  be  refcred  as  TP-TP  (Turning  Point  - 
Turning  Point). 

Let  us  express  the  solution  in  the  exponential  region  near  the  surface, 

0  ^  2<f,,  as 

Z„  = - -  sinh(07-(/>,) 

V  I  Kin  I 

where 


A 

4>T=\  \k! 


<^1=  I  \K2„\dz, 

Jq 

and  the  solution  in  the  exponential  region  near  the  bottom,  ij^z  <  H,  as 


■  cosh(</>2 -04,) 
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^b  =  \^  \K,n\dz- tanh  — 


_J _  d\'<!n\  \ 


In  addition,  let  us  express  the  solution  in  the  critical  region  centered  at  z  =  z,  as 


2  /  ^  \ 
^:n  ~  y\K^~^\)  * 


C\  ~  Y\ 


C,=-y!'Vz,). 

and  in  the  other  critical  region,  centered  at  z  =  Zj,  as 


K^„  =  y2(z2-^)  • 


r-  -1/6 

02  =  yin  yj 


By  letting 


C2  = 


16 


1 


Cl  =• 


-20, 


and  connecting  (23)  to  the  solution  in  the  oscillatory  region  i,<^2<^2j,  we  get 

1 


sin(<^t+-J-“i) 


(25) 


where 


and 


4>t  =  J 


a,  =  tan 


-1/  c 


2e^‘ 


Note  that,  in  obtaining  (25),  the  trigonometric  identity 


a  cos  d+b  sin  6  = -^1  a^+b^  sin^0+ tan 


has  been  used.  In  the  same  way,  by  letting 


(26) 


C2=' 


-2<5» 


and  connecting  (24)  to  the  oscillatory  region  solution  we  get 

1 


2.= 


sin(4>2  +  +1x2'^ 


where 


(27) 


<t> 


2  =  f  Vrn 


\dz 
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and 


Realizing  that  the  two  solutions.  (25)  and  (27),  in  the  oscillatoo’  region 
must  be  identical,  we  obtain 

=  sin(07+-J+a2)  .  (28) 

Using  (28),  the  characteristic  equation  for  this  class  of  normal  modes  is  obtained: 


I  Y^rt+a,-a2  . 

To  compute  the  normal  mode  we  can  use  (20)  for  z  ^  ij  and  (17)  for  z  ^  Zj. 
d.  Class  IV:  Pressure  Release  -  Rigid  Bottom 

Finally,  we  must  consider  modes  having  no  turning  points.  We  call  this  class 
PR'RB  (Pressure  Release  -  Rigid  Bottom).  Here  we  express  the  solution,  which  is  always 
oscillatory,  as 


with 


oc  — =r  sin  0+ 


cos  (}> 


(29) 


K^Jdz  . 


Since  the  surface  boundary  condition  must  be  satisfied,  we  must  have  i  =  0,  and  hence 


Z„  ~  — sin  0  . 


(30) 


On  the  other  hand,  (30)  must  satisfy  the  bottom  boundary  condition  (14)  and  this  leads 
to  the  characteristic  equation  equation  for  the  PR  -  RB  modes: 
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CH 

K^^dz  =  nn—  tan 

•'0 

for  Z)  >  0,  and 


Mi 

K^ndz  =  nn+  tan 

Jq 

for  Z)  <  0,  where  D  is  given  in  (22). 

2.  Formulae  Associated  With  Near-Boundary  Turning  Points 

In  our  derivation  of  the  formulae  in  the  previous  section,  we  have  applied  the 
boundary  conditions  to  the  oscillatory  or  exponential  region  solutions,  assuming  that  the 
boundaries  are  nowhere  close  to  any  turning  point.  In  the  special  case  that  a  bound.ary 
lies  inside  a  critical  region,  the  respective  boundary  condition  must  be  applied  to  the 
critical  region  solution  instead.  DilTerent  formulae  associated  with  the  first  three  classes 
of  modes  for  this  special  case  will  be  derived  next. 
a.  Class  I 

PR  -  TP  modes  have  lower  turning  points.  Since  the  solution  in  the  critical 
region  must  now  satisfy  the  bottom  boundary  condition  we  obtain 


where 


and 


At'  and  Bi'  are  the  derivatives  of  the  Airy  Functions  with  respect  to  C  and  are  defined 
by  (with  C  ^  0) 
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- -jr  ^  T «'”) 

AiV)  -  -  y  c['-!»(  y  ('")-'2n{  y  {''')'j 


5/'(a^ — Wc 

v'3 


/-.3(|C='’)+/.3(|C='^)' 


Defining 


and 


5,  =  /J/'(C„) 


and  applying  the  connection  formulae,  the  solution  in  the  oscillatory  region  becomes 


A 

y/K. 


wmt  /A  *  MBS  1 

1’’+  4  r 

^C0S[4>+-^J 

2n 


V  '^2n 


An  application  of  the  surface  boundary  condition  gives  the  following  characteristic 
equation; 


b.  Class  II 

TP  -  RB  modes  have  upper  turning  points.  To  satisfy  the  surface  boundary 
condition,  we  require  the  solution  in  the  critical  region  to  vanish  at  the  surface,  i.  e., 

Z„  =  fl/(C5V^/(C)-^/K5)<^5/(C) 


with 


1/3/' 


C5  =  y 


and 
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Defining 


A,  =  BiiHs) 


and 


and  applying  the  connection  formulae,  the  solution  in  the  oscillatory  region  becomes 
Z„  =  — ^sin(0+-^)-— 5^cos(0+-^)  . 

\  ^2n  V  ^2/1 

The  subsequent  application  of  the  bottom  boundary  condition  leads  to  the  following 
characteristic  equation  : 


for  D  <0,  and 


dz  — 


Ai-¥DBi 

B2-DA2 


A1ADB2 

B2-DA2 


for  D  >  0,  where  D  is  defined  in  (22). 
c.  Class  HI 

With  the  above  results,  we  can  easily  that  only  two  modifications  in  the 
general  formulae  for  TP  -  TP  modes  are  required.  These  include  replacing  by  /ij  and 

-y-  by  5j,  if  the  upper  turning  point  is  close  to  the  surface,  and  e**  by  /4,  and - — 

by  5„  if  the  lower  turning  point  is  near  the  bottom. 

3.  Criterion  For  Formulae  Selection 

We  must  define  a  criterion  for  switching  from  the  general  formulae  to  the  for¬ 
mulae  associated  with  near-boundary  turning  points.  It  was  found  by  Bender  and 
Orszag  (1978)  that  the  critical  region  extends  on  each  side  of  the  turning  point  until 
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Id  ~  1  .  In  accordance,  we  will  use  the  general  formulae  unless  the  distance  between  a 
turning  point  and  a  boundarj'  corresponds  to  values  of  |d  smaller  than  1. 

4.  Normalization 

Once  normal  mode  parameters  are  computed  using  the  appropriate  formulae 
including  the  characteristic  equations,  the  normalization  of  the  normal  modes  can  be 
achieved  by  numerical  integration  over  depth.  The  normalized  modes  (Z„'s)  are  related 
to  the  Z„  's  by 

z  =  c  z 


with 


1 


Zldz  . 

‘’o 


5.  Parameter  Storage 

We  now  define  the  necessary  parameters  required  to  completely  characterize 
each  mode.  The  first  parameter  is  obviously  the  class  number.  By  checking  the  formulae 
developed  in  the  previous  sections  it  is  seen  that  all  the  constants  {4>„  </>*,  D,  etc.)  can 
be  computed  from  the  horizontal  wavenumber,  the  depths  of  the  turning  points  and  the 
depth  of  the  ocean.  Therefore,  the  parameters  required  to  parameterize  each  class  of 
modes  are: 

a.  Classes  I  And  II 

1.  Class  number 

2.  Horizontal  wavenumber 

3.  Normalization  constant 

4.  Depth  of  the  turning  point 

b.  Class  III 

1 .  Class  number 

2.  Horizontal  wavenumber 

3.  Normalization  constant 

4.  Depth  of  the  upper  turning  point 

5.  Depth  of  the  lower  turning  point 


c.  Class  IV 


1.  Class  number 

2.  Horizontal  wavenumber 

3.  Normalization  constant 
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III.  ACCURACY  TEST 


A.  EXACT  ANALYTICAL  SOLUTIONS 

To  quantify  the  accuracy  of  the  first  order  WKB  approximation,  we  will  compare 
the  WKB  solutions  to  exact  analytical  solutions  to  the  Helmholtz  Equation  for  three 
abstract  sound  speed  profiles.  The  exact  solutions  are  derived  in  this  section.  We  use 
as  boundary  conditions  pressure  release  surface  and  rigid  bottom  for  all  three  cases. 

1.  Positive  Exponential  Profile 

In  this  upward  refractive  profile,  sound  speed  increases  exponentially  with  the 
depth  and  is  given  by 

c{z)  =  CQe^^ 

where  /?  is  a  constant.  The  surface  is  at  z  =  0  and  the  bottom  at  z  =  H.  The  Helmholtz 
Equation  is 


with  kJ  =  . 

^’o 

With  a  change  of  variable 


(31)  becomes 


After  division  by  and  using 


(32) 
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and 


a 


n 


(33) 


we  obtain 


This  is  the  Bessel  Equation  and  its  general  solution  is 

Z„  =  aJ^JiiXoX)+bY^S'^oX)  ■ 

The  boundary  conditions  are 

Z„(2  =  0)^Z„(x=l)  =  0 


dZ„  dZf. 


The  application  of  (35)  and  (36)  to  (34)  results  in  two  algebraic  equations 

«4;ao)+fey;>o)  =  o 


(34) 


(35) 

(36) 


(37) 


«-/',>o;t//)+^J''J«oXH)  =  0  .  (38) 

As  this  linear  system  of  algebraic  equations  is  homogeneous,  there  is  a  nontrivial  sol¬ 
ution  only  if  the  corresponding  Jacobian  is  zero,  i.e., 

•/»>o)  Y4<x,)J'4aaii)  =  0  . 

This  characteristic  equation  must  be  solved  in  order  to  obtain  the  a„'s,  which  are  the 
scaled  eigenvalues.  After  the  a„'s  are  evaluated,  the  k„'%  can  be  determined  using  (33).  It 
follow’s  from  the  surface  boundary  condition,  expressed  in  (37),  that  the  specific  solution, 
aside  from  a  multiplicative  constant,  is 

z,(2)oc  • 


The  constant  of  proportionality  is  given  by  normalization. 
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The  general  solution  to  (40)  is 

Using  the  same  anolog^'  as  in  the  previous  case,  we  obtain  as  characteristic 
equation 

=  0  (41) 

and  specific  solution 

Z„(z)oc  y,;ao)7,>o/")-7,>o)n„(“o«^0  • 

The  K„’s  can  be  found  using  (33)  after  solving  (41)  for  the  a,'s. 

3.  Hyperbolic  Profile 

In  this  third  case  we  assume  that  the  sound  speed  has  a  minimum  at  z  =  Zq.  The 
sound  speed  profile  is 

c(z')  =  Cq  cosh()?z') 
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with 


2  =  2—2 


The  Helmlioltz  Equation  to  be  solved  is  now 


dz 


,2 


•  + 


(t) 


2  1  2 

- 3 - 

cosh\Pz') 


Z„  =  0 


(42) 


With  the  following  change  of  coordinates 

X  =  tanhipz') 


(42)  becomes 


2  dZ„ 


+['<o(i-;(V'<yz„  =  o 


(43) 


After  dividing  (43)  by  )J^( !—;(’)  and  using  the  defmitions  given  in  (32)  and  (33),  (43)  can 
be  recast  as 


d 

dx 


il-x)' 


dx 


+  0^0- 


1-x' 


Zn-0 


(44) 


This  is  the  Associated  Legendre  Equation.  The  general  solution  to  (44)  is 

Z,  =  aP':{x)+bQ':{x)  , 


with 


-I  = 


and 


v(v+l)  =  ao  . 

The  functions  P“{x)  and  Qi(x)  are  the  Associated  Legendre  Functions  of  the  first 
and  second  kind  of  order  n  and  degree  v  .  Applying  the  boundary  conditions  we  get  the 
system  of  equations 


<‘iXs)+t>Q::iXs)  =  0 


(45) 
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a/’'v(;t//)+^<?':a;/)  =  0 


with 


Xs  =  tanh(-/?zo) 


and 


XH  =  tanh[)3(//-zo)]  . 

Thus,  the  corresponding  characteristic  equation  is 

p':ixs)Q'':(XH)-Q:iy.s)nixH) = o .  (46) 

Solving  (46)  for  fi  gives  the  eigenvalues  a„ 's  and  hence  the  horizontal  wavenumbers 
X,  's.  Using  (45)  the  solution  becomes 

Z„(2')  oc  Q:’-(x5}P:''(  tanh  Pz')-P:^{xs)Q:''i  tanh  fiz')  . 

B.  ANALYTICAL  EVALUATION  OF  WKB  PHASE  INTEGRALS  AND 
DERIVATIVES 

Now  we  use  the  WKB  formulae  to  obtain  the  first  order  WKB  solution  for  the  three 
analytical  sound  speed  profiles.  The  objective  in  this  section  is  to  evaluate  analytically 
all  the  related  WKB  integrals  and  derivatives  so  that  numerical  errors  can  be  eliminated 
in  one  of  the  error  analysis.  Horizontal  wavenumbers  computed  using  both  the  exact  and 
approximate  (i.e.,  WKB)  characteristic  equations  will  be  compared  in  the  next  section. 
The  evaluation  of  the  horizontal  wavenumbers  is  emphasized  because  they  are  the  key 
parameters  in  normal  mode  computations. 

1.  Positive  Exponential  Profile 

The  vertical  coordinate  z  is  positive  downward,  the  surface  is  at  z  =  0  ,  the  bot¬ 
tom  at  z  =  //  and  the  sound  speed  profile  is 

c(z)  =  Coe^^  . 


With 


the  corresponding  wavenumber  profile  is 
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K  =  K^e 


Let  us  first  consider  modes  that  have  a  turning  point  at 
(oscillatorv’  region),  is  given  by 


=  2.  ror0<2:<z 


with 

Kq 

The  spatial  phase  <f)  in  this  region  is  given  by 


=  J  \K2n\<i2 


With  the  change  of  coordinate 


(47)  becomes 


y.{2) 


In  the  region  z  ^2  ^  H  (exponential  region),  we  must  use 


krJ  =\j’xl-Kle 


The  spatial  phase  is  now 


(47) 


0 


=  \^\^zn\d2  , 


or,  with 
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+yn 


We  must  also  compute 


1 


2kinl 


2 


d\K,„\ 

dz 


which  can  be  equated  as 


Z)  = 


kIp 

2 


Figure  I  shows  a  PR  -  TP  mode  in  this  upward  refractive  medium. 
If  the  normal  mode  is  PR-RB  the  phase  is  given  by 


<f)  = 


or 


y„  tan 


Xiz) 

xiO) 


with 


r.U)  -  e-^‘ 


and 


4P 

2  ■ 

For  each  case  we  must  also  apply  the  respective  characteristic  equation. 

2.  Negative  Exponential  Profile 

For  this  profile  we  have  two  types  of  normal  modes:  TP  -  RB  and  PR  -  RB. 
The  sound  speed  is 
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Figure  1.  PR  -  TP  Mode 


c(z)  =  ro^  . 


Consequently,  for  a  TP  -  RB  normal  mode  and  for  z  ^  f  (oscillatory'  region)  we  have 

=  'f0\/  • 


So  in  this  region  the  phase  is 


and  D  is  given  by  (48).  Again,  for  each  of  the  classes  the  respective  characteristic 
equation  must  be  used. 

3.  Hyperbolic  Profile 

We  consider  two  types  of  normal  modes:  TP  -  TP  and  PR  -  RB.  For  the  type 
TP  -  TP  and  in  the  oscillatory  region  (i',  ^  2  <  z'j)  the  phase  will  be 
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and 


(49)  becomes 


0  =  «o 


For  the  exponential  region  near  the  bottom  (H>z'>  z'j)  the  phase  is 


4,= 


or 


an 

1+^ 

- 

2 

in 

^+sjU'-vl 

4- in 

VnJ -i^  +  O+yl 

- 

where 


For  this  type  of  normal  mode  D  is  given  by 

2[KWo{i-e„)r  ' 

In  the  exponential  region  near  the  surface  (z',  >  z')  the  phase  is 
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< 


«() 

In 

-l-y„  In 

1+^  y..v -n-<?)+y^ 

-1 

2 

- 

If  the  mode  is  PR  -  RB  the  phase  will  be 


0  =  oto 


sin 


-1 


I— sin  ’ 


Vn\J 


and  D  is  given  by 

Figure  3  shows  a  TP  -  TP  mode  trapped  by  the  refractive  index. 

C,  COMPARISON  OF  RESULTS 

1.  Horizontal  Wavenumber  Error  Analysis 

The  exact  characteristic  equations  derived  in  section  A  for  the  three  abstract 
profiles  were  solved  using  iterative  procedures  to  obtain  the  benchmark  k„ 's.  Similarly, 
the  approximate  WKB  characteristic  equations  developed  in  the  previous  chapter  were 
also  solved  for  the  three  profiles  with  the  use  of  the  algebraic  equations  developed  above. 
Very  low  frequencies  were  chosen  for  the  comparison  intentionally.  This  would  give  the 
WKB  method  a  real  test,  since  WKB  is  inherently  a  high-frequency  approximation.  For 
the  exponential  profiles,  we  used  <w  =  10  5"'  or  a  frequency  of  1.59  IIz.  For  the 
hyperbolic  profile  the  frequency  used  was  even  lower,  it  was  0.23  Hz.  Some  of  the 
computed  horizontal  w'avenumbers  (k„)  are  presented  in  Tables  1  through  5.  The  unit  for 
K„  is  inverse  meter  (w').  The  depth  of  the  ocean  is  taken  to  be  10000  m. 

In  Tables  1  and  2  the  exact  and  WKB  results,  as  well  as  the  absolute  and  relative 
errors,  for  the  positive  exponential  profile  are  displayed.  Tables  3  and  4  show  the  results 
for  the  negative  exponential  profile.  In  Table  5,  the  exact  and  approximate  results  for 
the  hyperbolic  profile  are  compared.  At  0.23  Hz,  mode  2  is  TP  -  TP  and  mode  3  is  PR 
-  RB. 

Overall,  we  can  see  that  the  absolute  error  (the  absolute  value  of  the  WKB  re¬ 
sult  minus  the  exact  result)  varies  between  10'^  and  10'’  m  '  and  the  relative  error  (the 
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Figure  3.  TP  -  TP  Mode 


absolute  error  divided  by  the  exact  value)  varies  between  10  *  and  10'*.  It  is  noted  that 
for  the  modes  with  turning  points  not  in  close  vicinity  to  the  boundaries  the  absolute 
error  usually  stays  below  10‘^  nr'.  The  fact  that  the  error  increases  when  a  turning 
point  is  close  to  a  boundary  is  easily  explained.  The  solutions  in  the  critical  regions  are 
the  less  accurate  because  they  use  a  linear  approximation  for  k^,.  So  when  they  are  used 
to  match  the  boundary  conditions  the  error  increases. 
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Table  1.  POSITIVE  EXPONENTIAL  PROFILE,  PR  -  TP  MODES 


Mode 

# 

Exact  K„  (w*) 

WKB  K„  (m-‘) 

Absolute 
Error  (/«-') 

Relative 

Error 

5 

.43861  X  10-' 

.43866  X  10-' 

5  X  10-’ 

l.l  X  10-- 

6 

.40441  X  10-' 

.40445  X  10-' 

4  X  10-’ 

1.0  X  lO-- 

7 

.37227  X  10-' 

.37230  X  10-' 

3  X  10-’ 

.8  X  10-- 

8 

.34179  X  10-' 

34182  X  10-' 

3  X  10-’ 

.9  X  10-- 

9 

.31274  X  10-' 

.31277  X  10-' 

3  X  10-’ 

1  X  10-- 

10 

.28556  X  10-' 

.28568  X  10-' 

1.2  X  10-‘ 

4.2  X  10-- 

Table  2.  POSITIVE  EXPONENTIAL  PROFILE,  PR  -  RB  MODES 


Exact  K„  (m-‘) 

WKB  K,  (m-') 

Absolute 
Error  (w-') 

Relative 

Error 

12 

.23190  X  10-' 

.23041  X  10-' 

1.5  X  10-* 

6.3  X  10-' 

13 

.18524  X  10-' 

.18490X  10-' 

1.8  X  10-' 

14 

.10756  X  10-' 

.10736  X  10-' 

2  X  10-‘ 

1.9  X  10-' 

Table  3.  NEGATIVE  EXPONENTIAL  PROFILE,  TP  -  RB  MODES 


Exact  K„  (m-‘) 

WKB  K„  (m-‘) 

Absolute 
Error  (m-‘) 

Relative 

Error 

23 

.82307  X  10-' 

.82308  X  10-' 

1  X  10-’ 

1  X  10-' 

24 

.79463  X  10-' 

.79464  X  10-' 

I  X  10-’ 

1  X  10-' 

25 

.76664  X  10-' 

.76665  X  10-' 

1  X  10-’ 

1  X  10-' 

26 

.73903  X  10-' 

.73904  X  10-' 

1  X  10-’ 

1  X  10-' 

27 

.71158  X  10-' 

.71155  X  10-' 

3  X  10-’ 

4  X  10-' 

In  general,  for  each  type  of  mode  the  WKB  approximation  gets  better  as  mode 
number  increases.  The  only  exception  is  when  turning  points  are  close  to  the  boundaries. 
To  see  this,  let  us  examine  the  results  for  the  positive  exponential  profile  (Tables  1  and 
2).  Starting  from  the  lowest  modes  that  have  one  turning  point,  the  error  decreases  as 
mode  number  increases.  At  mode  number  10,  the  accuracy  of  the  WKB  method  de- 
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Table  4.  NEGATIVE  EXPONENTIAL  PROFILE,  PR  -  RB  MODES 


Mode 

# 

Exact  K„  (nr') 

WKB  >c„  (nr') 

Absolute 
Error  (m-') 

Relative 

Error 

29 

.65250  X  10-^ 

.65091  X  10-2 

1.6  X  10-* 

2.4  X  10-2 

30 

.61746  X  10-^ 

.61683  X  10-2 

1.6  X  10-5 

1.0  X  10-2 

31 

.57723  X  10-2 

.57688  X  10-2 

3.5  X  10-‘ 

6.1  X  lO--* 

32 

.53082  X  10-2 

.53058  X  10-2 

4.5  X  10-* 

33 

.47670  X  10-2 

.47652  X  10-2 

1.8  X  10-‘ 

3.8  X  10-* 

Table  5.  HYPERBOLIC  PROFILE 


Exact  K„  (nr') 

WKB  K„  (m-‘) 

Absolute 
Error  (nt-') 

Relative 

Error 

2 

.901 14  X  10-2 

.89881  X  10-2 

2.3  X  10-‘ 

2.2  X  10-2 

3 

.74139  X  10-2 

,72953  X  10-2 

1.2  X  10-2 

1.6  X  10-2 

creases  because  the  turning  point  comes  too  close  to  the  bottom  boundary.  For  mode 
12  and  higher  a  turning  point  does  not  exist  and  the  oscillatory  region  is  bounded  by  the 
physical  boundaries.  After  the  change  of  mode  class  the  error  starts  to  decrease  again 
as  mode  number  increases. 

Errors  in  the  k„'$  for  all  cases  are  small  enough  that  they  do  not  significantly 
influence  the  vertically  integrated  phases  0  because  the  maximum  ocean  depth  is  of  the 
order  of  km's.  In  the  horizontal  direction  the  horizontal  phase  is  approximately  (this 
phase  is  exact  when  c=  c(z)).  The  error  in  k„  limits  the  range  for  which  the  use  of  WKB 
is  accurate.  If  we  want  to  keep  the  phase  error  below  a  few  degrees,  the  tolerable  error 
in  K„  is  of  the  order  of  10"*  m-‘  for  a  range  of  10  km,  10'‘  m-'  for  100  km  and  lO"’  w' 
for  1000  km.  All  the  k„'s  associated  with  the  modes  that  do  not  have  interactions  with 
the  bottom  boundary,  as  calculated  from  the  WKB  method,  have  errors  less  than 
10~*  m-‘,  implying  that  the  method  is  good  to  at  least  100  km  at  1  Hz.  For  higher  fre¬ 
quencies  the  results  would  be  much  better. 

2.  Errors  In  The  Interference  Distances 

The  acoustic  pressure  square  at  far  field  can  be  WTitten  as  (assuming  no  range 
dependence) 
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where  the  A„'s  are  constants  and  the  acoustic  source  is  at  z  =  Zg  and  r  =  0  (Clay  and 
Medwin,  1977).  We  can  see  that  the  interference  terms  are  functions  of  the  differences 
of  horizontal  wavenumbers.  The  interference  distance  (or  interference  w'avelength)  de¬ 
fined  by 
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nm 
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is  therefore  an  important  parameter  for  transmission  loss  calculation  (Chiu  and  Ehret, 
1990).  In  general,  dominant  interferences  are  between  adjacent  modes  (Chiu  and  Ehret, 
1990).  So  let  us  see  what  is  the  size  of  the  errors  in  Ak„  =  k„|  computed  using  the 

WKB  approximation.  In  Tables  6  and  7  we  display  the  exact  and  WKB  Ak:,'s,  as  well 
as  the  absolute  errors.  As  expected,  we  can  see  that  the  errors  in  Ak:,  var>'  in  a  similar 
W'ay  as  in  k„  and  are  slightly  smaller. 


Table  6.  POSITIVE  EXPONENTIAL  PROFILE,  PR  -  TP  MODES 


n 

Exact  Ak„  (wi*‘) 

WKB  Ak„  (w-*) 

Absolute  Error 
(/«-') 

5 

.3420  X  10-5 

.3421  X  10-3 

1  X  10-’ 

6 

.3214  X  iO-5 

.3215  X  10-3 

1  X  10-' 

7 

.3048  X  10-’ 

.3048  X  10-3 

0 

8 

.2905  X  10-3 

.2905  X  10-3 

0 

9 

.2718  X  10-3 

.2709  X  10-3 

9  X  10-’ 

3.  General  Numerical  Method 

The  W’KB  results  in  the  previous  section  were  obtained  using  analytical  means 
rather  than  numerical  methods.  The  reason  for  that  was  we  wanted  to  see  how  accurate 
the  WKB  approximation  is  regardless  of  the  numerical  methods  used  to  evaluate  inte¬ 
grals  and  derivatives.  As  we  see,  absolute  errors  in  the  order  of  10-’  m"'  for  the  k„'s  are 
achieved  at  1  Hz.  This  means  that  the  WKB  solution  is  very  accurate,  especially  because 
we  used  ver>'  low  frequencies. 
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Table  7.  POSITIVE  EXPONENTIAL  PROFILE,  PR  -  RB  MODES 


n 

Exact  Ak„  (m-‘) 

WKB  Ax,  («r>) 

Absolute  Error 
(m-') 

12 

.4666  X  10-^ 

.4551  X  10-5 

1.15  X  10-s 

13 

.7768  X  10-5 

.7754  X  10-5 

1.4  X  10-‘ 

For  an  arbitrary  profile  all  the  integrations  and  diflerentiations  must  be  done 
numerically.  It  is  expected  that  the  errors  will  increase  due  to  numerical  noise.  The 
Fortran  program  WKBGEN  (Appendix  C)  developed  in  this  thesis  can  be  applied  to  an 
arbitrary  profile.  For  each  mode  it  finds  class,  horizontal  wavenumber  and  depths  of  the 
turning  points.  This  program  was  also  applied  to  the  three  abstract  profiles.  Some  nu¬ 
merical  results  for  the  ic„'s  are  presented  in  Tables  8  through  11. 

Tables  8  and  9  show  the  exact  and  numerical  WKB  results  and  the  respective 
absolute  error:  for  the  positive  exponential  profile.  Tables  10  and  11  show  to  the  cor¬ 
responding  rc:”>'.a  for  the  negative  exponential  profile. 


Table  8.  POSITIVE  EXPONENTIAL  PROFILE,  PR  -  TP  MODES 


Exact  K„  (w-*) 

WKB  X,  (w-') 

Absolute 
Error  (m-') 

5 

.43861  X  10-5 

.43874  X  10-5 

1.3  X  10-‘ 

6 

.40441  X  10-5 

.40454  X  10-5 

1.3  X  10-‘ 

7 

.37227  X  10-5 

.37234  X  10-5 

7x  10-’ 

8 

.34179  X  10-5 

.34184  X  10-5 

5  X  10-’ 

9 

.31274  X  10-5 

.31284  X  10-5 

1  X  10-‘ 

10 

.28556  X  10-5 

.28574  X  10-5 

1.8  X  lO-" 

As  we  can  see,  the  errors  in  the  k/s  computed  using  the  numerical  method  are  slightly 
larger  than  the  previous  results  (Tables  1  through  4)  but  are  approximately  in  the  same 
order  of  magnitude. 
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Table  9.  POSITIVE  EXPONENTIAL  PROFILE,  PR  -  RB  MODES 


Mode 

u 

Exact  <„  (//i  ’) 

WKB  K,,  (m-') 

Absolute 
Error  (m-') 

12 

.23190  X  10-^ 

.23044  X  10-2 

1.5  X  10-* 

13 

.18524  X  10-2 

.18494  X  10-2 

3.0  X  10-* 

14 

.10756  X  10-2 

.10744  X  10-2 

1.2  X  10-‘ 

Table  10.  NEGATIVE  EXPONENTIAL  PROFILE,  TP  -  RB  MODES 


Mode 

a 

Exact  K„  (nr') 

WKB  K„  (m~') 

Absolute 
Error  (/«-') 

23 

.82307  X  10-2 

.82309  X  10-2 

2  X  10-’ 

24 

.79463  X  10-2 

.79469  X  10-2 

6  X  10-’ 

25 

.76664  X  10-2 

.76669  X  10-2 

5  X  10-^ 

26 

.73903  X  10-2 

.73909  X  10-2 

6  X  10-’ 

27 

.71158  X  10-2 

.71159  X  10-2 

1  X  10-' 

Table  11.  NEGATIVE  EXPONENTIAL  PROFILE,  PR  -  RB  MODES 


Mode 

Exact  K,  (»i-') 

WKB  K„  (m-‘) 

Absolute 
Error  (w-') 

30 

.61746  X  10-2 

.61689  X  10-2 

5.7  X  10-‘ 

31 

.57723  X  10-2 

.57689  X  10-2 

3.4  X  10-‘ 

32 

.53082  X  10-2 

.53059  X  10-2 

2.3  X  10-‘ 

33 

.47670  X  10-2 

.47649  X  10-2 

2.1  X  10-‘ 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  WKB  approximations  of  acoustic  normal  modes  seem  to  give  results  accurate 
enough  to  be  of  practical  use.  Although  WKB  is  inherently  a  high  frequency  approxi¬ 
mation,  meaning  that  it  works  better  for  higher  frequency,  the  results  from  our  tests 
show  that,  even  for  frequencies  around  one  Hertz,  this  technique  has  an  appreciable 
accuracy. 

When  applying  the  WKB  algorithm  for  arbitrary  sound  speed  profiles  the  determi¬ 
nation  of  the  class  of  each  normal  mode  must  be  done  carefully.  Each  class  has  different 
formulae.  There  are  a  total  of  four  classes  in  a  single  duct  or  channel  environment. 

A  small  w’eakness  of  the  method  is  that  the  error  in  the  WKB  solution  increases 
when  a  turning  point  is  very  close  to  a  boundary.  However,  the  corresponding  errors  are 
still  very  small  for  the  exponential  and  hyperbolic  profiles  used  in  this  study. 

The  exact  analytical  solutions  to  the  Helmholtz  Equation  can  also  be  used  for 
comparison  to  any  other  methods.  These  exact  solutions  are  subject  to  pressure  release 
surface  and  rigid  bottom  boundary  conditions.  Exact  solutions  can  also  be  found  for 
any  boundary  conditions.  Specifically,  we  can  maintain  the  pressure  release  surface 
boundary  condition,  which  is  a  good  assumption,  and  use  a  non-rigid  bottom  boundary 
condition. 

Some  difficulties  were  found  when  working  with  Bessel  and  Associated  Legendre 
Functions.  The  IMSL  subroutines  used  for  Bessel  Function  evaluation  cannot  handle 
some  orders  and  values  of  the  argument.  With  respect  to  the  evaluation  of  Associated 
Legendre  Functions,  subroutines  are  not  available  in  the  IMSL  libraries.  Some  Fortran 
programs  were  coded  to  compute  them.  These  programs  also  work  only  for  certain 
ranges  of  order,  degree  and  arguments  of  the  functions.  Further  programming  work 
concerning  these  transcendental  functions  is  recommended. 

WKB  formulae  for  modes  trapped  in  the  water  column  over  a  non-rigid  bottom  were 
derived  in  Appendix  B  although  they  are  not  implemented  for  this  analysis.  A  general 
mathematical  expression  for  boundary  conditions  for  normal  modes  is  also  presented  in 
Appendix  B.  The  use  of  this  mathematical  expression  is  not  a  trivial  matter  because  the 
expression  depends  on  the  density  and  the  sound  speed  profile  of  the  sediment.  Further 
studies  on  using  WKB  algorithms  for  arbitrary  bottom  boundary  conditions  are  recom¬ 
mended. 
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In  this  thesis  only  single  duct/channel  environments  have  been  considered.  Essen¬ 
tially,  we  have  ignored  double  duct/channel  problems.  However,  the  WKB  approxi¬ 
mation  method  has  been  applied  to  other  nonacoustic  but  equivalent  double 
duct/channel  problems  (for  example,  transmission  of  electromagnetic  waves  through 
potential  barriers)  with  sucess  (Bender  and  Orszag,  1978). 

In  conclusion,  the  WKB  method  allows  for  accurate  and  fast  computation  of  normal 
modes  and  their  eigenvalues.  In  addition,  it  allows  for  storage  of  the  results  in  terms  of 
only  a  few  parameters  for  each  mode. 
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APPENDIX  A.  FIRST  ORDER  WKB  THEORY 


With 


cl4>  =  \K^„\dz 


the  Helmholtz  Equation  becomes 

d^Z„  1  d\Kj  dZ„ 

d(^^  |k^„|  dij> 


(A.l) 


This  equation  has  some  similarities  with  a  Bessel  Equation.  Let  us  try  a  solution  in  the 
form 


Z„  oc 


Substituting  this  trial  solution  form  in  (A.l)  we  get 


d^F  1  dF 
d<p^  4> 


l+r(2)- 


i\l2? 


F«0 


(^.2) 


with 


riz)=l-h 


^zn  I 
dz 


1  1 

2  dz^ 


If  r(r)  is  negligible,  i.e.,  |  r{z)  1  <11,  (A.2)  can  be  approximated  as 


d^F  1  dF 

d<f>'^  <t>  d<i> 


1- 


0' 


f=0 


w’hich  is  a  Bessel  Equation. 

The  general  solution  to  (A.3),  for  >  0,  is 

F(</>)  =  a'y,p(<^)+fe'r,/2(0) 


(<4.3) 


{AA) 


where  J,,2  and  y,,^  are  the  order  1/2  Bessel  Functions  of  the  1"  and  2"''  kind,  respectively. 
But  since 
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(A.4)  becomes 


and  thus 


,  X  sin  0  .  ,  cos  0 
f(0)  =  a— ^+6 

V  0  V 


^  ^  <2  sin  0+^)  cos  0 

Z„(z)  = - == - 


If  kI,  <  0,  instead  of  (A. 3)  we  would  get  the  Modified  Bessel  Equation  whose  general 
solution  is  given  by 

/•(0)  =  a'ii:,/2(0)+6'/,/2(0) 


where  A'l/j  and  /i,j  are  the  Modified  Bessel  Functions.  Since 


we  now  have 


A^,/2(0)oc-^=r 

V  0 


/,/2(0)  cx: 


F(0)  = 


and  thus 


Z„{2)  = 


ae  *+be'^ 
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Let  us  assume  that  there  is  a  turning  point  at  z  =  0,  k^>0  for  r>0,  for 

<  0  and  near  the  turning  point  kI,  =  yz.  So,  close  to  the  turning  point  we  have 


4>  =  j  \K,„\dz  =  s-jy'''^\z\ 


3/2 


with  5  =  z/ 1 2 1 ,  and  after  some  algebra, 

1  1 


Using  this  in  (A.l)  we  get 


\kJ  d<p  30 


- -  +  ~ - —+Z  =0 

^  30  d<f> 


Let  us  now  look  for  a  solution  near  the  turning  point  in  the  form 

Z„oc(^^<t>y^F{4>). 

Substituting  this  solution  form  in  (A. 5)  we  get 


dLL  + 

1 

JL. 

d<t>^ 

<t> 

dit> 

1- 


(1/3) 

</>' 


2  n 


f=o 


(^.5) 


whose  general  solution  is 


f(0)  =  C,5_,/3(0)+C25i/3(0) 

w’here  B  represents  Bessel  Functions.  For  xf,  >  0  they  arc  and  Ji,y  For  <  0  they 
are  /.,/3  and  /,/j.  Relating  Z„  -with  F  we  obtain 

Z,  =  (y  0)'''Cc,B_,/3(0)+C25,/3(0)] 


or 


(.4.6) 
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where  the  plus  sign  is  for  kJ,  <  0  and  the  minus  sign  for  kI„  >  0.  A  way  to  avoid  the  in¬ 
convenience  of  sign  switching  is  to  use  Airy  Functions  which  aie  related  to  5^.,  3  (see 
Chapter  III,  Section  A). 
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APPENDIX  B.  GENERAL  BOUNDARY  CONDITIONS 
A.  GENERAL  FORMULA 

Until  now  we  have  assumed  that  the  boundary  conditions  are  pressure  release  at  the 
surface  and  rigid  bottom.  The  first  one  is  a  good  assumption  but  the  second  one  can  be 
far  away  from  reality.  Let  us  now  present  a  general  bottom  condition  formula  for 
normal  modes  that  depends  on  the  bottom  characteristics. 

The  acoustic  pressure  due  to  the  mode  is 

p„  =  Z„iz)R„{r,d)e“^‘  .  (5.1) 

The  vertical  component  of  the  momentum  equation  relating  p„  to  the  vertical  particle 
velocity  in  the  n"'  mode  is 


Pi 


Vn 


Oi 


SPn 

CZ 


(5.2) 


where  />,  is  the  water  density  and  the  bottom  is  assumed  to  be  fiat.  Using  (B.l)  in  (B.2) 
we  get 


dZ„ 


lu}{ 


(5.3) 


At  the  bottom  boundary,  the  requirement  is  that  both  p„  and  Uy„  are  continuous 
across  the  water- sediment  interface.  In  other  words,  the  normal  acoustic  impedance 
across  the  interface  is  continuous.  The  normal  acoustic  impedance  associated  with  the 
n*  mode  is 


Pn_ 

Sn 


(5.4) 


By  using  (B.l)  and  (B.3)  in  (B.4)  one  can  obtain  as  a  general  condition 

f  dZ„  _  ,  ^  \ 


dz 


+<■ 


•iVn 


=  0 


(5.5) 


if  Zf,„  at  the  water- sediment  interface  is  known  from  the  sediment  properties. 
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Let  us  use  a  local  plane  wave  approximation  of  the  normal  mode  near  z  —  II  in  the 
water  column,  i.e.. 

Z„{z)  = 

where  is  a  constant,  R  is  the  complex  plane  wave  reflection  coefTicient  and  k„„  is  the 
vertical  wavenumber  in  water  computed  at  the  boundary’  (Clay  and  Medwin,  1977). 
With  this  plane  wave  approximatio,  one  can  easily  show  that  (Kinsler  et  al.,  1982) 


where  the  index  1  refers  to  water,  2  refers  to  sediment  and  b  refers  to  the  value  at  the 
boundary.  R  is  also  known  as  the  Rayleigh  Reflection  Coefficient  which  can  be  equated 
as 
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and,  for  a  pressure  release  boundary-  (i?  =  —  1), 

(Z,W  =  0  . 

B.  TRAPPED  MODES  IN  THE  WATER  COLUMN 

We  will  consider  only  modes  that  are  trapped  in  the  water  layer  in  this  WKB  for¬ 
mulation.  This  means  that  is  purely  imaginary  or 

Klzn  =  iPn 


with  P„  real  and  defined  by 


(5.9) 


Substituting  (B.9)  in  (B.7)  we  get 


dz 


fL 

Pi 


PnZn 


=  0 


(5.10) 


Since  (B.7)  is  a  general  boundary  condition,  (B.IO)  is  also  general  but  is  only  appli¬ 
cable  to  modes  trapped  in  the  water  column.  We  will  apply  (B.IO)  to  each  class  of 
normal  modes  in  the  water  column. 

1.  Class  I 

As  usual  we  assume  that  there  is  a  turning  point  at  z  =  f  with  0  <  z  <  //.  In  the 
exponential  region  (z  <  z  <  //)  the  solution  is 

- J=cosh(«^-</>t)  (5.11) 

V 


with 


and 


^zn  I 
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I  I  / .  2  W 

I  ^2n  I  2 


As  (B.ll)  must  satisfy  (B.IO),  ^4  must  be 


f".  .  .  .-1/  I  .Pi  Pn  \ 

^»=J.  2Ik:.I=  ■*  l-<„l  /”« 


(B.12) 


with  /?,  given  by  {B.9).  The  only  difference  relative  to  the  rigid  bottom  case  is  expression 
for  4>b-  All  the  other  formulae  remain  unchanged. 

2.  Class  II 

In  the  oscillatory  region,  i  <  2  <  //,  the  solution  is 

f )  (S.13) 

\/  K,-,  ^?>i 


=  J  \K,„\dz 


<^*5=  I  ' 

Jq 


(B.13)  must  satisfy  (B.IO),  implying  that 


r  r 

sinU;;+-j  +  tan  '  - 

I  l-L^-iD+E}e^- 


(5.14) 


where  D  is  defined  in  Chapter  II,  Equation  (22), 
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\2=li  • 


(B.15) 


£  = 


P\  Pn 
Pi 


Therefore,  the  characteristic  equation  following  from  (B.14)  is 


dz  = 


TC—OC 


for  D<0,  and 


for  D>0,  where 


dz  = 


n—a 


a 


tan 


.n±§.g-*,  -1 

-(D+£)e^'  - 


All  the  other  formulae  are  identical  to  the  rigid  bottom  case. 

3.  Class  III 

For  this  class  the  only  difference  relative  to  the  rigid  bottom  case  is  that  0*  must 
now  be  computed  with  (B.12). 

4.  Class  IV 

The  solution  over  the  entire  ocean  depth  is 

Z„  =  — ==-  sin  <t> 
v'Xzn 


with 


To  satisfy  (B.IO)  we  require 


5? 


dz 


=  nn—  tan 


for  D>0,  and 


dz  =  nn+  tan 


for  £><  0,  with  E  given  by  (B.15).  This  is  now  the  Class  IV  mode  characteristic  equation. 
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APPENDIX  C.  FORTRAN  PROGRAM  VVKBGEN 


A.  PROGRAM  DESCRIPTION 

The  program  WKBGEN  is  formed  by  a  main  program  and  several  subprograms. 
The  main  program  evaluates  the  minimum  sound  speed  and  then  controls  the  subpro¬ 
grams.  The  inputs  must  be  given  in  the  PARAMETER  statement.  The  inputs  include 
the  ocean  depth  (H),  acoustic  frequency  (F),  the  wavenumber  increment  used  in  the  k„ 
search  (DK),  first  mode  (NMI)  and  last  one  (N.MF)  to  be  computed.  The  units  in  this 
program  are  those  of  the  MKS  system.  The  input  sound  speed  profile  is  specified  using 
the  subprogram  FUNCTION  C(Z).  SUBROUTINE  TYPE  evaluates  the  class  of  each 
mode.  SUBROUTINE  CHARAC  controls  the  search  for  k„  and  the  turning  points.  The 
appropriate  characteristic  equation  is  given  by  FUNCTIOxN  EQCHAR  and  the  turning 
points  are  evaluated  by  SUBROUTINE  ZTURN.  The  auxiliary  subprograms  FUNC¬ 
TION  PHASE,  FUNCTION  FKZ2  and  FUNCTION  FKZ  evaluate  respectively  the 
phase  integral,  and  |k„|.  The  Airy  Functions  are  computed  by  the  FUNCTION'S 
AI(Z)  and  BI(Z)  and  their  derivatives  by  FUNCTION'S  DAI(Z)  and  DBI(Z).  These  four 
subprograms  use  the  SUBROUTINE'S  MMBSJR  and  MMBSIR  in  the  IMSL  libraries 
to  evaluate  the  Bessel  and  Modified  Bessel  Functions,  respectively. 

As  output,  the  mode  numbers  and  the  respective  eigenvalues  are  printed  on  the 
screen.  These  results,  as  well  as  the  depths  of  the  turning  points,  are  also  written  in  a  file 
whose  name  is  specified  in  the  CALL  EXCMS  statement  at  the  begining  of  the  main 
program. 

The  numerical  method  used  to  solve  for  the  characteristic  equation  and  to  find  the 
turning  depths  is  the  simple  but  safe  Bisection  Method  (Gerald  and  Wheatley,  1989). 
The  method  to  compute  integrals  is  the  Trapezoidal  Rule  (Gerald  and  Wheatley,  1989). 
Derivatives  are  evaluated  by  forward  or  backward  finite  differences  (Gerald  and 
Wheatley,  1989). 

B. 

C 
C 

c* 
c 
c 
c 
c 


PROGRAM  LISTING 
PROGRAM  WKBGEN 


MAIN  PROGRAM 
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o  cj  n  non  ono  o  o  o  o  nn 


IMPLICIT  REAL*8  (A-H,0-Z) 

PARAMETER  (H=  ,F=  ,DK=  ,NMI=  ,NMF=  ) 


CALL  EXCMSC'FILEDEF  1  DISK  WKBGEN  DATA  A') 

0M=8.  D0*DATAN( 1.  D0)*F 

ZEX=DM0D(H,10.D0) 


CS=C(0) 

CB=C(H) 


CM=99999.D0 
DO  Z=0.D0,H-ZEX,10.D0 
CM=DMIN1(CM,C(Z)) 
END  DO 

CM=DMIN1(CM,C(H)) 


XKO=OM/CM 

XKL=XK0 


DO  N=NMI,NMF 
C 

DO  XKI=XKL,0,DK 
C 

XKF=XKI+DK 

CNI=OM/XKI 

CNF=OM/X'KF 

CALL  TYPE(CNI,CS,CB,NTYPE) 

CALL  ZTURN (CNI,H,ZEX,NTYPE,ZI1.ZI2) 

CALL  ZTURN( CNF , H , ZEX , NTYPE , ZFl , ZF2 ) 

CALL  CHARAC ( OH , H , ZEX , N , XKI , XKF , NTYPE ,ZI1,ZI2,ZF1,ZF2,NS0L,XKS0L, 
$  ZT1,ZT2) 

IF  (NSOL.EQ.  1)  GO  TO  100 
C 

END  DO 
C 

100  PRINT*, N.XKSOL 
C 

IF  (NTYPE.  EQ.  1.  OR.  NTYPE.  EQ.  2)  THEN 
WRITE  (1,*)  N, NTYPE, XKSOL.ZTl 
ELSE  IF  (NTYPE.  EQ.  3)  THEN 

WRITE  (1,*)  N, NTYPE, XKS0L,ZT1,ZT2 
ELSE 

WRITE  (1,*)  N, NTYPE, XKSOL 
END  IF 
C 

XKL=XKS0L 


55 


o  o  o  o  o  o  o  o  o  n  n  o  n  n  n  n 


END  DO 

STOP 

END 


Jr**'*********ifc***'*-************************************'**'***'***********iV* 


SUBROUTINE  TYPE(CN,CS,CB,NTYPE) 
REAL*8  CN,CS,CB 

IF  (CN.GT.CS.  AND.CN.lt.  CB)  THEN 
NTYPE=1 

ELSE  IF  (CN.  LT.CS.AND.cn.  GT.CB)  THEN 
NTYPE=2 

ELSE  IF  (CN.  LT.CS.AND.cn.  LT.CB)  THEN 
NTYPE=3 
ELSE 

NTYPE=4 
END  IF 

RETURN 

END 


SUBROUTINE  CHARAC ( OM , H , ZEX , NM , XI , XF , NT , Z I 1 , Z I 2 , ZF 1 , ZF2 , NSOL , XSOL , 
$  Z21,Z22) 

IMPLICIT  REAL*8  (A-H,0-Z) 

XSOL=9999. DO 
NSOL=l 

FI=EQCHAR(NM,NT,H,OM,XI,ZIl,ZI2) 

FF=EQCHAR( NM , NT , H , OM , XF , ZF 1 , ZF2 ) 

IF  ((FI*FF).GT.  O.DO)  THEN 
NSOL=0 
GO  TO  100 

ELSE  IF  (FI.  EQ.  O.DO)  THEN 
XSOL=XI 
GO  TO  500 

ELSE  IF  (FF.EQ.  O.DO)  THEN 
XSOL=XF 
GO  TO  500 
ELSE 

CONTINUE 
END  IF 

DO  200  N=l,500 
XM2=(XF+XI)/2.D0 
IF(N. EQ. 1)  GO  TO  50 
IF  ((XM2-XM1).EQ.  0. DO)  THEN 
XS0L=XM2 
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GO  TO  500 
END  IF 
50  XM1=XM2 
CN2=0M/XM2 

CALL  ZTURN( CN2 , H , ZEX , NT , Z2 1 , Z22 ) 
F2=EQCHAR( NM , NT , H , OM , X2 , Z21 , Z22) 
IF  ((F2*FF).  LT.  0.  DO)  THEN 
XI=XM2 
ELSE 

XF=XM2 
FF=F2 
END  IF 
200  CONTINUE 
100  CONTINUE 
500  CONTINUE 
RETURN 
END 


SUBROUTINE  ZTURN( CN , H , ZEX , NT , ZTl , ZT2 ) 
IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*8  NM 

ZT1=99999.D0 

ZT2=99999.D0 

DO  100  D=  0.D0,H,10.D0 

XF=D+10.  DO 

XI=D 

FI=C(XI)-CN 

FF=C(XF)-CN 

IF  ((FI*FF).GT,  O.DO)  THEN 
GO  TO  100 

ELSE  IF  (FI.  EQ.  O.DO)  THEN 
XSOL=XI 
GO  TO  500 

ELSE  IF  (FF. EQ. O.DO)  THEN 
XSOL=XF 
GO  TO  500 
ELSE 

CONTINUE 
END  IF 

DO  200  N=l,500 
XM2=(XF+XI)/2.D0 
IF(N.  EQ.  1)  GO  TO  50 
IF  ((XM2-XM1).EQ. 0. DO)  THEN 
XS0L=XM2 
GO  TO  500 
END  IF 
50  XM1=XM2 

F2=C(XM2)-CN 

IF  ((F2*FF).LT.  0.  DO)  THEN 
XI=XM2 
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ELSE 

XF=XM2 
FF=F2 
END  IF 
200  CONTINUE 
100  CONTINUE 
500  CONTINUE 
ZT1=XS0L 

IF  (NT.EQ.  3)  THEN 

DO  110  D=  H,0.D0,-10.D0 

XF=D-10.D0 

XI=D 

FI=C(XI)-CN 

FF=C(XF)-CN 

IF  ((FI*FF),GT. O.DO)  THEN 
GO  TO  110 

ELSE  IF  (FI.  EQ. O.DO)  THEN 
XSOL=XI 
GO  TO  510 

ELSE  IF  (FF.EQ.  O.DO)  THEN 
XSOL=XF 
GO  TO  510 
ELSE 

CONTINUE 
END  IF 

DO  210  N=l,510 
XN2=(XF+XI)/2.D0 
IF(N.EQ.  1)  GO  TO  51 
IF  ((XM2-XM1).EQ.  O.DO)  THEN 
XS0L=XM2 
GO  TO  510 
END  IF 
51  XM1=XM2 

F2=C(XN2)-CN 

IF  ((F2*FF).LT.  O.DO)  THEN 
XI=XM2 
ELSE 


XF=XM2 
FF=F2 
END  IF 
210  CONTINUE 
110  CONTINUE 
510  CONTINUE 
ZT2=XS0L 
END  IF 
RETURN 
END 
C 
C 
C 


C 

C 

C 


FUNCTION  EQCHARC NM , NT , H , OM , XI , ZTl , ZT2 ) 
IMPLICIT  REAL*8  (A-H,0-Z) 
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PI=4.D0*DATAN(1.D0) 

D=(FKZ(OM,XI,H)-FKZ(OM,XI,H-5.DO))/(10.DO*FKZ(OM,XI,H)**2) 
IF  (DABS(D).GT.  .999D0)  THEN 
D1=(D/DABS(D))*.  999D0 
ELSE 
D1=D 
END  IF 


IF  (NT.EQ. 1)  THEN 

EQCHAR=PHASE(O.DO,ZT1,OM,XI)-(NM-. 25D0)*PI 
IF  (H-ZTl.GT.  10.  DO)  THEN 
PHIB=PHASE ( ZT 1 , H , OM , XI ) -DATANHC D 1) 
EQCHAR=EQCHAR+DATAN(DEXP( -2.  D0*PHIB)/2.  DO) 

ELSE 

GA.M=(FKZ2(0M,XI,H)-FKZ2(0M,XI,ZT1))/(2T1-H) 
ZH=(GAM/DABS(GAM))*DABS(GAM)**(1.  D0/3.D0)*(H-ZT1) 
A=DBI(ZH) 

B=DAICZH) 

EQCHAR=EQCHAR-DATAN( B/A) 

END  IF 


ELSE  IF  (NT.  EQ.  2)  THEN 

EQCHAR=PHASE ( ZT 1 , H , OM , XI ) 

IF(ZT1.GT.  lO.DO)  THEN 
PHIS=PHASE(O.DO,ZT1,OM,XI) 

ALF 1=DEXP (PHIS) +D*DEXP( - 1 .  DO*PHI S ) / 2 .  DO 
ALF2=((DEXP(-1.D0*PHIS)/2.D0)-D*DEXP(PHIS)) 

IF  (DLOG10(DABS(ALF1))-DLOG10(DABS(ALF2)).GT. 60.D0)  THEN 
IF  (D.  GE.O)  THEN 
ALF=PI/2.  DO 
ELSE 

ALF=-1.D0*PI/2.D0 
END  IF 
ELSE 

ALF=DATAN( ALF 1 / ALF2 ) 

END  IF 

IF  (D.  GE.O.  DO)  THEN 

EQCHAR=EQCHAR-(NM-1. 25D0)*PI+ALF 
ELSE 

EQCHAR=EQCHAR-(NM-.  25D0)*PI+ALF 
END  IF 
ELSE 

GAM=(FKZ2(OM,XI,ZT1)-FKZ2(OM,XI,O.DO))/ZT1 
ZS=(GAM/DABS(GAM) )*DABS(GAM)**( 1.  DO/3.  DO) 

A=BI(ZS) 

B=AI(ZS) 

IF  (D.  LT.  O.DO)  THEN 

EQCHAR=EQCHAR- ( NM- .  25D0 )*PI+DATAN( ( A+D*B ) / ( B -D*A) ) 

ELSE 

EQCHAR=EQCHAR- ( NM- 1 .  25D0 )*PI+DATAN( ( A+D*B ) / ( B -D*A) ) 

END  IF 
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c 

c 

c 


END  IF 


C 

C 

C 


C 

C 

C 


ELSE  IF  (NT.  EQ. 3)  THEN 
IF  (H-ZT2.GT.  lO.DO)  THEN 
PHIB=PHASE(ZT2,H,OM,XI)- DATANH( D 1 ) 
ALF2=DATAN(DEXP( -2.  D0*PHIB)/2. DO) 

ELSE 

GAM=(FKZ2(0M,XI,H)-FKZ2(0M,XI,ZT2))/(ZT2-H) 
ZH=(GAM/DABS(GAM))*DABS(GAM)**(1. DO/3.DO)*(H-ZT2) 
A=DBI(ZH) 

B=-l.DO*DAI(ZH) 

ALF2=DATAN(B/A) 

END  IF 

IF(ZT1.GT.  lO.DO)  THEN 
PHIS=PHASE( 0. DO , ZTl , OM ,XI ) 

ALF1=DATAN(DEXP( -2.  D0*PHIS)/2.  DO) 

ELSE 

GAM=( FK22 ( OM , XI , ZTl ) -FKZ2( OM , XI , 0.  DO ) ) / ZTl 
ZS=(GAM/DABS(GAM) )*DABS(GAM)**( 1.  DO/3.  D0)*ZT1 
A=BI(ZS) 

B=AI(ZS) 

ALF1=DATAN(B/A) 

END  IF 

EQCHAR=PHASE( ZTl , ZT2 ,OM,XI) 

EQCHAR=EQCHAR-(NM-.  5D0)*PI-ALF1+ALF2 


ELSE 

EQCHAR=PHASE ( 0 .  DO , H , OM , XI ) 

IF  (DABS(D).GT.  l.D-60)  THEN 
E=l.  DO/D 

IF(D.  GT.  0.  DO)  THEN 

EQCHAR=EQCHAR-NM*PI+DATAN( E ) 
ELSE 

EQCHAR=EQCHAR-NM*PI -DATAN( E ) 
END  IF 
ELSE 

EQCHAR=EQCHAR-NM*PI+PI/2.  DO 
END  IF 


END  IF 


RETURN 

END 


C 

C 

C 

C*********************************************************************** 


c 

c 

c 
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FUNCTION  PHASE(A,B,OM,XKN) 

IMPLICIT  REAL*8  (A-H,0-Z) 
XNL=DINT((B-A)/5.D0) 

NL=NINT(XNL) 

EX=DMOD((B-A),5.DO) 

PHASE=0.  DO 
DO  1=1, NL 
XI=DBLE(I) 

PHASE=PHASE+(FKZ(OM,XKN,(XI-l.DO)*5.DO+A) 
$  +FKZ(OM,XKN,XI*5.DO+A))*2. 5D0 

END  DO 

PHASE=PHASE+(FKZ(OM,XKN,XNL*5.  DO+A) 

$  +FKZ(0M,XKN,B))*EX/2.D0 

RETURN 
END 


C 

C 

C 


C 

C 

C 

FUNCTION  FKZ(OM,XKN,Z) 

REAL*8  FKZ,OM,XKN,Z 
FKZ=(OM/C(Z))**2'XKN**2 
FKZ=DSQRT( DABS( FKZ ) ) 

RETURN 

END 

C 

C 

c 

c*********************************************************************** 


c 

c 

c 

FUNCTION  FKZ2(OM,XKN,Z) 

REAL*8  FKZ2,OM,XKN,Z 

FKZ2=(OM/C(Z))**2-XKN**2 

RETURN 

END 


C 

C 

C 


C 

C 

C 

FUNCTION  C(Z) 

REAL*8  C,Z 
C= 

RETURN 

END 

C 

c 

c 

Q********'************************'********  ******************************* 
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FUNCTION  AI(Z) 

IMPLICIT  REAL*8  (A-H,0-Z) 
REAL*8  J13,JM13,I13,IM13 
DIMENSION  RJ(2),WK(4),B(2) 


IF  (Z.EQ.  O.DO)  THEN 
AI=. 35502805D0 
RETURN 

ELSE  IF  (Z.  LT.  O.DO)  THEN 


ARG=2.  D0*(DSQRT( -1.  D0*Z)**3)/3.  DO 
N=2 


ORDER=2.DO/3.  DO 

CALL  MMBS JR( ARG , ORDER , N , R J , WK , lER ) 
JM 13=4.  DO*R J( 1 ) / ( 3 .  DO*ARG ) -RJ( 2 ) 


ORDER=1.DO/3.DO 

CALL  MMBS JR( ARG , ORDER , N , R J , WK , lER ) 
J13=RJ(1) 


AI=DSQRT( -1. DO*Z)*( JM13+J13)/3.  DO 


RETURN 

ELSE 


ARG=2. D0*(DSQRT(Z)**3)/3.  DO 

ORDER=2.  DO/3.  DO 

NB=2 

IOPT=l 


CALL  MMBS IR( ARG , ORDER , NB , lOPT , B , lER) 
IM13=4. DO*B( l)/(3. D0*ARG)+B(2) 


ORDER=1.DO/3.DO 
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CALL  MMBSIRC ARC, ORDER, NB,IOPT,B,IER) 
I13=B(1) 


AI=DSQRT(Z)*(IM13-I13)/3.D0 

RETURN 


END  IF 

END 


fr********************************************************************** 


FUNCTION  BI(Z) 

IMPLICIT  REAL*8  (A-H,0-Z) 
REAL*8  J13,JM13,I13,IM13 
DIMENSION  RJ(2),WK(4),B(2) 


IF  (Z.EQ.  O.DO)  THEN 
BI=. 61492663D0 
RETURN 

ELSE  IF  (Z.  LT.  0.  DO)  THEN 


ARG=2.  D0*(DSQRT( -1.  D0*Z)**3)/3.  DO 
N=2 


ORDER=2.DO/3.DO 

CALL  MMBSJRC ARC, ORDER, N.RJ.WK.IER) 
JM13=4. D0*RJ( l)/(3. D0*ARG)-RJ(2) 


ORDER=1.DO/3.DO 

CALL  MMBS JR( ARC , ORDER , N ,RJ , WK , lER) 
J13=RJ(1) 


BI=DSQRT( -1. D0*Z/3. D0)*( JM13-J13) 
C 
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RETURN 

ELSE 


ARG=2. D0*(DSQRT(Z)**3)/3.  DO 

ORDER=2.DO/3.DO 

NB=2 

I0PT=1 


CALL  MMBS IR( ARC , ORDER , NB , lOPT , B , lER) 
IM13=4. D0*B( l)/( 3. D0*ARG)+B( 2) 


ORDER=1.DO/3.DO 

CALL  MMBS IR( ARG , ORDER , NB , lOPT , B , lER) 
I13=B(1) 


BI=DSQRT(Z/3.D0)*(IM13+I13) 

RETURN 


END  IF 


END 


FUNCTION  DAI(Z) 

IMPLICIT  REAL*8  CA-H,0-Z) 
REAL*8  J23,JM23,I23,IM23 
DIMENSION  RJ(2),WK(4),B(2) 


IF  (Z.  EQ.  O.DO)  THEN 
DAI=-. 25881940D0 
RETURN 

ELSE  IF  (Z.  LT.  O.DO)  THEN 
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ARG=2. D0*(DSQRT( -1.  D0*Z)**3)/3.  DO 
N=2 


0RDER=1. DO/ 3.  DO 

CALL  MMBSJR( ARG, ORDER, N,RJ,WK,IER) 
JM23=2. D0*RJ( l)/(3.  D0*ARG)-RJ(2) 


ORDER=2.DO/3.  DO 

CALL  MMB  S JR( ARG , ORDER , N , R J , WK , lER) 
J23=RJ(1) 


DAI=Z*(JM23-J23)/3. DO 


RETURN 

ELSE 


ARG=2. D0*(DSQRT(Z)**3)/3.  DO 

ORDER=1.DO/3.DO 

NB=2 

I0PT=1 


CALL  MMBS IR( ARG , ORDER , NB , lOPT , B , lER) 
IM23=2.  D0*B( l)/(3.  D0*ARG)+B(2) 


0RDER=2.  D0/3.D0 

CALL  MMBS IR( ARG , ORDER , NB , lOPT , B , lER) 
I23=B(1) 


DAI=-1. D0*Z*( IM23-I23)/3.  DO 
RETURN 


END  IF 


END 
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C 

C*********  *******************  ***-fr***'<r*  A' vie  ***il(*******ifr'i<n>r'Jc'Jr'Jr'i»r****^r^']lr'******'ft 


FUNCTION  DBI(Z) 

IMPLICIT  REAL*8  (A-H,0-Z) 
REAL*8  J23.JM23,I23,IM23 
DIMENSION  RJ(2),WK(4),B(2) 


IF  (Z.EQ.  O.DO)  THEN 
DBI=. 44828836D0 
RETURN 

ELSE  IF  (Z.LT.  O.DO)  THEN 


ARG=2.  D0*(DSQRT( >1.  D0*Z)**3)/3. DO 
N=2 


ORDER=1.DO/3.DO 

CALL  MMBSJRCARG, ORDER, N,RJ,WK,IER) 
JM2  3=2 . D0*R J( 1 ) / ( 3 . D0*ARG ) -R J( 2 ) 


0RDER=2.  DO/3.  DO 

CALL  MMBS JR( ARG , ORDER , N , R J , WK , lER) 
J23=RJ(1) 


DBI=- 1. D0*Z*( JM23+J23) /DSQRT( 3.  DO ) 


RETURN 

ELSE 


ARG=2 . D0*( DSQRT( Z ) **3 ) / 3 .  DO 

ORDER=1.DO/3.DO 

NB=2 

I0PT=1 


CALL  MMBS IR( ARG , ORDER , NB , lOPT , B , lER ) 
IM23=2. D0*B( 1 )/( 3. D0*ARG)+B( 2) 
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0RDER=2. DO/ 3. DO 

CALL  MMBSIRC ARC, ORDER, NB.IOPT.B.IER) 
I23=B(1) 


DBI=Z*(IM23+I23)/DSQRT(3.D0) 

RETURN 


END  IF 


END 
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